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form factor, skin friction, heat transfer, and velocity profiles. 

This report includes a listing of the FORTRAN IV computer program with an explana- 
tion of the equations involved, the solution methods, and the preparation of input. Numer- 
ical examples are included to illustrate the input and show what type of results are ob- 
tained. Running times are less than 1 minute on IBM- 7094 equipment. 
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is given of the input required and the output received. Numerical examples are included 
to illustrate the input and show what type of results are obtained. 

The report is divided into two main parts: INFORMATION FOR GENERAL USER 
and INFORMATION FOR PROGRAMMER . Those wishing to use the program in its 



present form need only read the first part, INFORMATION FOR GENERAL USER . This 
contains descriptions of the general method, the preparation of input, and the output re- 
ceived. Information of interest to someone who may wish to change the program is con- 
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describing the individual subroutines. 




laminar methods have. It is one of the most accurate, programmable, general methods 
available for the laminar case. 

(6) The Sasman-Cresi method (ref. 5) is used for the solution of the turbulent bound- 
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The following statements briefly describe some of the input and output features of the 
program: 
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The following sections discuss the differential equations solved in the laminar and 
turbulent sections of the program. The methods used in predicting transition and sepa- 
ration are also outlined. 



After the program reads input data and does preliminary calculations, the laminar 



= N(n, 


transformed free-stream velocity, 





For isothermal, or nearly isothermal, surfaces, the solution of equation (1) is simplified 
since N can be expressed as a linear function of the pressure gradient parameter (cor- 
relation number) n as follows: 



When appropriate initial values of n are used, equation (4) is solved numerically along 
the boundary-layer surface, giving the distribution of pressure gradient parameter (or 
correlation number ) n. 

Once n is obtained at each station, the other boundary-layer and heat-transfer 
parameters are easily obtained. The shear-stress parameter l and the heat-transfer 




equation is obtained by applying a Mager-type transformation (ref. 13) to Frandtl’s equa- 
tions in which flow variables appear as time- averaged quantities. The momentum equa- 



and 




transformed free-stream velocity, 
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Equations (8) and (9) are the coupled, first-order, ordinary differential equations 
which govern the development of the turbulent boundary layer. These equations are not 
uncoupled in the Sasman-Cresci method as they were in Reshotko and Tucker’s analysis. 

Using initial values calculated in the laminar routine or given by the user, equations 
(8) and (9) are solved by a Runge-Kutta scheme, giving the distribution of f and Hj 
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an NACA airfoil which has both laminar and turbulent boundary layers, and the second is 
a converging-diverging channel with turbulent boundary layer throughout. The airfoil il- 
lustrates the prediction of transition by the program, and the channel illustrates the use 
of initial values of displacement and momentum thickness at the beginning of the boundary 
layer. In both cases the computer input is included. 
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Figure 1. - Comparison of analytical results with experimental data on NACA 0012 airfoi 
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Figure 3 shows the input variables as they are punched on the data cards. The first 
input data card is for a title. In any of the columns from 1 to 80 of this card, the user 
may put whatever information he wishes in order to identify the data deck. The remain- 
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Dictionary of Input Variables 


After a title card, the following input variables are given: 
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KATCH integer (0 or 1) indicating whether laminar-boundary-layer separation (if 
encountered) should reattach as a turbulent boundary layer: 

Separation and stop 

Reattach 



real variable used when KATCH = 1, indicating ratio of momentum thickness 
after reattachment to momentum thickness at laminar separation 



of velocity and Mach number is desir 


integer (0 or 1, seeKPRE) indicating whether printing of solutions of lami- 
nar and turbulent differential equations is desired (see OUTPUT) 
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array of free-stream Mach numbers relative to surface atX Y input 
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Surface geometry . - The X and Y input arrays give coordinates of points which de- 
scribe the surface geometry (see figs. 4 and 5). These are the points at which pressure 
(velocity, etc. ) and wall temperature are given, and at which boundary-layer parameters 
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nation point at station 1 (KLE = 0), or a sharp or pointed leading edge at station 1 
(KLE = 1). 

The program will also accept data which begin at some point within the boundary 
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11 me user gives a value to NTURB, he may or may not choose to give values to 
DTURB and TTURB. If he does so, they are used at the station specified by NTURB. If 
values are not given to DTURB and TTURB, laminar values of displacement and momen- 
tum thickness at the station specified by NTURB are used as initial values. (When 
NTURB = 1, however, values must be given for DTURB and TTURB. ) 



Use of initial values of 5* and 6. - The user may give initial values whenever he 
knows them. In the laminar case, only one value is required (either DLAM or TLAM). 
If both are given, TLAM is used. If initial turbulent values are given, both DTURB and 
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Mass, kg (or N-sec /m) 

Temperature, K 

With these consistent sets of units, density is mass per unit volume, energy is force 
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Instability and transition stations are located no matter how far RTHI has gone past 
RCRIT or RTRAN, respectively. Examination of the values of RCRIT and RTRAN some- 
times indicates that instability or transition occurs closer to the previous station than the 
indicated station. Since final answers for turbulent thickness parameters are affected by 



the station at which transition occurs, the user may wish to rerun the program with addi- 
tional points in this area so as to locate transition more exactly. 
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(2) THERE IS A STAGNATION POINT AT A STATION OTHER THAN STATION 1. 
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This message is printed by subroutine LAMNAR if NTURB = 1 and values are not given 
for both DTURB and TTURB. If the turbulent boundary layer is to start at station 1, ini- 
tial values must be given for both displacement thickness and momentum thickness. 
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This message is printed by subroutine LAMNAR if the value of correlation number calcu- 
lated in the solution of the laminar-boundary-layer equation is beyond the range -0.32 < 
n < 0. 16. This can be caused by a strong adverse or favorable pressure gradient, gen- 
erally occurring rapidly over one or two stations along the surface. Several input sta- 
tions should be used to pass through any rapid change in pressure or velocity distribution. 


(10) IF INITIAL TURBULENT VALUES ARE GIVEN, 
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tines. The subroutine LAMNAR solves the laminar differential equation, computes the 
laminar-boundary-layer parameters, and checks for transition. After transition, it com 
putes the initial values for the turbulent boundary layer. TURBLN solves the turbulent- 
boundary-layer equations by using RUNKUT, and computes the turbulent-boundary-layer 
parameters. Finally, PROFIL computes both laminar and turbulent velocity profiles. 
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The following subroutines are particular to the BLAYER 



COMMON blocks described in the previous section. More general subroutines (i. e. , they 
are independent of BLAYER) are discussed in the section AUXILIARY SUBROUTINES. 
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puted from UE, and all five distributions are reprinted. After surface distributions have 
been finalized, other constant parameters depending on them are calculated. 

Calculation of surface gradients . - Surface gradients of velocity and Mach number 
(du g /dx, dM e /dx, and dMg/dfc/L)) are calculated by means of calls on either GRADNT 
or SPLINE, depending on which the user has requested. 
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These functional relations are given in appendix B. 

Initial values . - A check is made for initial values given by the user. If initial val- 
ues are given for the turbulent boundary layer at station 1, the program transfers to the 
end of LAMNAR, where transformed quantities are computed for use by the TURBLN 
routine. If an initial value is given for thelaminar boundary layer, the corresponding 



correlation number nj is computed, since the laminar equation is solved in terms of n 
If a laminar momentum thickness is given, the corresponding correlation number n is 
computed from the following relation obtained by rearranging equation (15), which is de- 
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ROOT, depending on whether the Mach number gradient dM e /d(x/L) is positive or nega 
tive. Since 5* = 6H and 6 and H are both functions of n, 8* can be expressed as a 
function of n. Subroutine ROOT solves the following equation for n 1 , given 6*: 
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A spacing of l/500th of surface length L is used to step from point to point. This 
spacing represents a compromise between use of machine time and excessive accuracy in 
the solution of the equation. The resulting surface lengths and correlation numbers are 
stored in tables. 



Calculation of laminar-boundary-layer parameters . Through interpolation, correla- 
tion numbers are obtained at the input-output stations on the boundary-layer surface. 
Using these values of n, the boundary-layer parameters are calculated. 

Momentum thickness is obtained from the following equation (see appendix A): 





Boundary-layer thickness 6 is obtained from equation (41) of reference 4: 
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and heat transfer per unit area is obtained from 
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The form factor is obtained from the incompressible form factor by the equation 
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local Nusselt number is computed as defined: 



Subroutine PROFIL prints all the principal bound ary -layer parameters computed by 



and the shape factor 
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which is part of equation (14). 
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The following general subroutines are called by the main subroutines, but do not use 
the COMMON variables. Each has its own variables, which are defined below instead of 
in the main dictionary. 


Subroutine SMOOTH is a simple data-smoothing routine. It smoothes an array 
length n by the following formula 


< 4-1 

O 

>> 


00 



in 









Figure 9. - Parameters used in calculation of gradients 
at input points. 


O) o 
a a 
O 

® S 

0 O 
& Ctf 


m + 

+ r- 


4- + 

•H 


+3 3 

m & 


00 

in 



Slopes at the end points are given by 


o> 

in 


00 


Oi 


o 

LO 


in 




i 






+ 

a 


ll 

$ 

o 

u 

a 


ca 

o 

ctf 


II 


B| 


Assuming a circle 


CM 

LO 


d 

0 

X 


CM 

03 


03 


cd 

n 

ca 


0 

O, 

o 

03 


0 

ha 

d 

cd 

X 

O 


X 

+ 

cd 


cd 


X 

+ 

cd 


0 

ft 

o 

i— H 

CQ 

d 

•rH 

0 

ha 

§ 

XI 

O 


03 

£ 

O 


03 
O S 

0 

fn 

Cd 

Q 

5 

O 

u 

o 


03 

0 

X 

cd 

•rH 

> 

H-> 

a 

d 

•rH 

0 

X 

H 


0 

3 

cd 

•rH 

u 

cd 

> 

H-H 

0 

T3 

0 

& 

73 

d 


XI 


B 

& 

fn 

fn 

cd 


03 

H-> 

g 

s 

0 

03 


fn 

0 

X 


d 

d 


XI 

o 


o 

0 

ft 


X 

X ft 


g 

•rH 

■§ 

fn 

ba 


B 

ft 

Q 


03 

£ 

O 


03 

cd 

0 

fn 

cd 

Eh 

£ 

Q 


fn 

S O 

£ fn 

^ o 

XI «*H 

03 

0 

° | 


fH 

cd 

> 

*3 

d 

fn 

0 


03 

oT 

d 

■rH 

O 

ft 

h-> 

g 

o 

cd 

•J— ■» 

TS 

cd 

d 

0 

0 


0 

X 

03 

0 

ft 

O 

03 


X 

GO 


03 


O 

ft 

<+-> 

g 

0 
cd 

•r— a 
■§ 

g 

0 

1 

X 

CQ 

0 

0 

g 

H-> 

CQ 


Eh 

a 

o 


D- 

OO 

a> 

c 

H— » 

o 


CO 


I 

cd 


cd 

d 

o 

ha 

cd 

•rH 

73 

• rH 

fn 


i 

>> 


73 

d 

o 

o 

0 

03 


O 0 
ft X 


X 

g 


0 

0 03 


cd 

X 

H-> 

d 

o 


r< -i-t H-> 

±> I> •«-* 


03 

0 

> 


I 

8 


73 

d 

o 

0 


03 o 
+n X 


ft X5 
d 
0 


0 
> 
fH 

1:3 d 

o S 


0 


0 


d .... 

•i-H «H 
r-H «4H 

ft 0 
03 O 


X 

d 

o 

0 

x 


0 S> 


03 

5 


w 

5 

X 

ft 

co 


© "Eo 

d 


§ 

fn 

X 

d 


d 

0 

1 

o< 


0 

X 


CQ 

03 0 

03 
d 


ft 

5 

x 

ft 

w CO 
d 
*d 

cd 0 


X 

o 


> 

fH 

d 

o 


C— *f— ! 


d 

• rH 

o 

ft 

0 

d 

•rH 

I*"*H 

ft 

03 

H— > 

d 

0 

o 

cd 

t? 

cd 

0 03 

f I 

rt -g 

H-> 64-H 

X $ 
h-> cd 

O 0 
H-» fn 

•a * 

d ft 

S' 5 


o 
d 
0 
fn 

03 0 o 

i-h <+h h; 

0 
5h 


03 


ft 

CQ 


ba 

d 

•rH 

> 

•rH 

ba 

d 

o 

0 

d 

.d 


co 55 .cd 

.a a 


ft 

a 

0 

fn 

0 

X 

•rH 

0 

H-> 

cd 

0 

> 

• rH 
H-> 

cd 

t> 

• rH 

fH 

0 

73 


X 

ft 

CO 

fn 

o 

<4-4 

03 

0 

X 

cd 

• rH 

fH 

cd 

s> 

H-> 

a 

d 

•rH 

0 

X 

ft 


CQ 

0 

H-> 

cd 

d 

•rH 

73 

fn 

O 

0 
0 

1 

X 


>> 

cd 

fn 

fH 

cd 


X 


ha 

d 

•rH 

73 

d 

o 

ft 

03 

0 

fn 

fn 

O 

0 

03 

0 

d 

Id 

> 

d 

o 

♦rH 

H-i 

0 

d 

d 


£ 

fH 

fH 

cd 


cd 

fH 

fH 

cd 


fH 

o 

X 


CQ 

H-> 

g 

s 

0 

'0 


fH 

0 

X 

s 

d 

d 


CQ 

£ 

o 


CQ 

cd 

0 

fH 

cd 

X 


X 

ft 

CO 


03 

0 

3 

cd 

•rH 

fH 

cd 
> 
H— > 

a 

H-> 

d 

o 

0 

x 

Eh 


03 

"d 

•rH 

o 

ft 

0 

d 

♦rH 

r-H 

ft 

03 

H— > 

cd 

X 


0 

0 

ft 

03 

0 

fn 


X >h Jz; 


«fH 

o 

03 

0 

► 

X 

cd 

> 

•rH 

fH 

0 

73 

H-> 

03 

fH 


>> 

cd 

fH 

fH 

cd 


X 

Q 

>H 

Q 


B 

d 

•rH 

o 

ft 

0 

d 

•rH 

ft 

03 

H— > 

cd 

X 


0 

0 

ft 

03 

0 

fn 


03 

0 

> 

•rH 

td 

> 

•rH 

fH 

0 

73 

"d 

O 

0 

0 

03 


>» 

cd 

fH 

fH 

cd 


CM 

B 

CM 

Q 


O 

O 

OcC 

CD 

c 

O 

-Q 

Z3 

on 


Subroutine ROOT locates a root for f(x) = y in the interval a, b (see fig. 11). The 
values of f(x) are calculated by another subroutine FUNCT. ROOT divides the interval 




Y y (see fig. 11) 

FUNCT external subroutine to calculate f(x) 
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value of independent variable at which value of dependent variable is desired 
The output variable for LGRNGE is 

value of dependent variable corresponding to ARG 



Subroutine CURVFT evaluates the polynomial f(x,y) using the coefficients C. The 


CM 






H 

T — 1 

1—1 

1 

1 

S 3 

1 

c 

c 

X 

X^ 

X 

S 3 


S 3 

CM 

CO 

o 

a 

U 

4 - 

+ 

4 - 



CM 



>> 

CM 


CM 

X 


X 

CO 

CM 

CO 

+ 

X 

4* 

S3 

CO 

S3 

CM 

U 

u 

u 

4- 

4- 

+ 




This subroutine is also used for evaluating functions of one independent variable. 
The input variables for CURVFT are as follows: 
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30 CONTINUE 

TEST ALL INTERVALS ASAIN IF ANY WERE SUBDIVIDED THE LAST T IHE 
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Using equation (A3), this becomes 
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Substituting equation (All) into equation (A9) yields 



DERIVATION OF EQUATION (20) 



(A13) 



DERIVATION OF EQUATIONS (18) AND (28) 


In reference 18 the following expression is derived, relating form factor H 
equivalent incompressible form factor H, for flows over insulated surfaces with 




it can be shown that 



Substituting equation (A19) into equation (A15) gives 
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CPR is used for the storage of these coefficients. 



(3) Nondimensional thermal conductivity of air as a function of nondimensional static 
temperature, (k/k .) = f(T/T 7 ): 



In LAMNAR the following curve fits are used: 

(1) Stagnation point correlation number against wall temperature function 
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curves in figure 3 of reference 4. The double- valued regions at the ends of the S 

W 

-0. 8, and -1.0 curves are not reflected in equation (B8). 

(6) Thickness ratio parameter against correlation number and wall temperature 
function , 6 tv/ 6 tr = f ( n > S J : 
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